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The Global Positioning System enables positioning in 3 dimensions about our planet. It has 
been operationnal since 1994. Twenty-four satellites are used to achieve this performance. The 
signals sent by these satellites are electromagnetic waves travelling through our atmosphere 
down to the small receivers used by the civilian community and the military. 

Because of varying meteorological conditions (namely temperature and humidity changes 
along the ray path), the rays do not travel in a straight line. They bend towards the surface. 
As a consequence, the ray path between two points is longer than a straight line, and the time 

it takes for a signal to travel this distance is longer. 

In 1995, a small GPS receiver was launched on a satellite (GPS/MET). It became possible 
to perform radio occultations around the Earth: the source - one of the 24 GPS satellites - is 
seen by the receiver as it rises or sets around the other side of the Earth. When the source 
disappears, the receiver progressively loses the signals. By measuring accurately the time delay 
between the emission and the reception of the signal, it is possible to infer which part of the 
delay is due to the atmosphere. 

We use GPS/MET data to retrieve temperature and humidity prohles simultaneously. 
A specific method is implemented : it combines information from numerical forecasts and 
GPS observations in an optimal way. Comparing the result with an independent source of 
observations (weather balloons), we demonstrate that GPS data have the potential to improve 
weather analyses. We also show that improved temperature and humidity profiles can be 
obtained using information from a forecast model. This confirms results obtained in this study 
using simulated data. 
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Abstract. 

The Global Positioning System (GPS) employs a constellation of satellites which provide a stable 
source of electromagnetic signals available for radio occultation purposes about our planet. The 
atmospheric-induced bending of the rays can be converted into refractivity profiles for each occultation. 
We implement here a one-dimensional variational (IDVAR) analysis which enables retrieving both 
temperature and humidity profiles without ambiguity. 

We perform a linear error analysis study and a fully non-linear Monte-Carlo simulation. They both 
show the potential impact of the GPS to improve upon an accurate first guess from a data assimilation 
system. The information content of the GPS soundings is significant at all latitudes for temperature, 
and only in the tropics for humidity. 

IDVAR retrievals based on real data collected during the GPS/MET experiment in 1995 are 
validated against radiosondes. They show an excellent capacity of the GPS measurements to resolve 
the tropopause. 

In the northern hemisphere and in the tropics, we demonstrate a net reduction of temperature 
bias and standard deviation, as compared with the Goddard Earth Observing System (GEOS) model 
background information used to constrain the profiles. The results show that an analysis using GPS 
refractivity can yield significant improvement over the background at a relatively small computing 
cost, and in spite of the spherical symmetry assumption used to derive refractivity from the GPS 
measurements. We also show that a IDVAR yields better results than a direct retrieval method in 
which a temperature (humidity) profile is assumed to be the truth in order to retrieve a humidity 
(temperature) profile. Finally, we point out significant sensitivities of the retrievals to the gravitational 
constant, the surface pressure and the number of levels used to perform the analysis. 

This represents a step towards using the GPS data to improve NWP forecasts. 


3 


1. Introduction 

The Global Positioning System (GPS) system enables positioning in 3 dimensions on our planet. 

Its radio signals are sensitive to the atmosphere. This makes it possible to perform soundings using the 
radio occultation technique (RO). This approach has been used at NASA for more than thirty years to 
study the atmospheres of other planets [e.g. Lmdal. 1992], Gorbunov and Sokolovskiy [1993] provided 
simulations of GPS radio occultation measurements. Kursinski et al. [1997] also simulated many 
aspects of GPS measurements and their expected error characteristics. GPS/MET (1995) was the first 
radio occultation experiment conducted on Earth using that technique [e.g. Kursmskz et al, 1996]. 

Before using the GPS data in a data assimilation system for numerical weather prediction (NVVP) 
or climate study, a first step is to appraise the impact these data have on the analysis. One technique 
for assessing impact is the Observing System Simulation Experiment (OSSE). Kuo et al [1998] have 
performed OSSE’s with simulated GPS data. There are known limitations of OSSE s [Atlas, 199 (] and 
therefore care must be taken when interpreting the results. For example, no observation errors were 
included in the Kuo et al [1998] OSSE. We will give an example of the effect of neglecting errors in a 

simulation. 

In this work, both simulated and real observations are analyzed to gauge the impact of GPS data 
on an analysis. For this purpose, we have developed a one dimensional variational (IDVAR) analysis. 
In this approach, background information from a General Circulation Model (GCM) forecast is used 
to constrain the retrievals. Analyzed temperature and humidity are obtained without ambiguity. The 
approach also yields an estimate of the errors, based on the assumed observation and background errors. 

The outline of the paper is as follows : First we give a description of the GPS. Next, the radio 
occultation technique is described. We discuss the different possible approaches to retrieve and 
assimilate atmospheric properties from GPS measurements. Then, we describe the implementation of 
a IDVAR analysis of refractivity. After this, we show results of a Monte-Carlo simulation and linear 
error analysis. We then discuss some unexpected sensitivities of the GPS analyses. Finally, we validate 
GPS/MET IDVAR retrievals using collocated radiosondes. 
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2. The Global Positioning System 

The GPS was designed by the U.S. Department of Defense (USDoD). Its main purpose is to aid m 
navigation. It has enabled accurate positioning anywhere about the globe since 1994. It consists of 24 
orbiting satellites distributed in six orbital planes. The 20,200 km orbits are circular with an inclination 
of 55® and a period of 12 hours. The principle of the GPS is based on the path delay of a propagating 
radio signal. Two L-band frequencies are used : Ly 1575.42 MHz (or 19.0 cm) and L. 1227.60 MHz (or 

24.4 cm). 

The accuracy of GPS measurements is intentionally degraded by encryption for security reasons. 
This is known as “selective availability”. First-generation receivers were launched on GPS/MET (1995) 
and Oersted (1999). Their ability to track the occultations is degraded when the encryption is on, 
which is generally the case except during very short periods for scientific studies. Second-generation 
GPS receivers built for the purposes of the radio occultation are to be flown on the upcoming missions 
Champ (Germany) and SAC-C (Argentina). These receivers are able to better process encoded signals, 

for which the Ln signal-to-noise (SNR) ratio is very low. 

An exact positioning cannot be achieved for several reasons [Doerflinger, 1998]. First, the exact 
position and velocity of the GPS satellites must be known perfectly (at a given time) to position a 
receiver relative to them. It follows that (1) orbit errors (2) clock drifts (3) relativistic effects (4) 
receiver system errors (antenna and receiver noise) and (5) use of approximate co-ordinates contribute 
to the total error of the measurement. Other perturbations in signals that degrade the quality of the 
positioning for a terrestrial receiver are caused by atmospheric refraction and scattering which may 
cause multipath effects. 

The applications of GPS to meteorology include positioning m situ measurements, computing 
wind speed from the drift of radiosondes, and measuring the total column of water vapor above 
a ground-based receiver [e.g. Bems et al., 1992, Elgered et ai, 1997]. In this study, we focus on 
space-based GPS receivers. 
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3, Radio Occultation with GPS 

3.1, Assets of the GPS Radio Occultation Technique 

There are several positive attributes of GPS occultations as an atmospheric sounding device. A 
high-inclination LEO provides a set of observations that covers the globe fairly uniformly at a relatively 
low cost [Yunck et al, 1988]. The GPS coverage anywhere about the globe is advantageous when 
comparing to that obtained from balloon-launched radiosondes (about 800 each 12 hours, the majority 
over the northern hemisphere continents). A single GPS receiver may optimally obtain approximately 
500 occultations per day [Kursmskt et al., 1997]. 

When compared within IR spaceborne sounders, the radio occultation (RO) technique with GPS 
has the advantage of being an ■ all-weather” system. It is relatively insensitive to aerosols, cloud or 
rain, due to wavelengths of order 20 cm. Unlike other techniques the GPS radio occultation provides 
a degree of self-calibration because relative phase shifts are the relevant information. Moreover, the 

stability of the transmitters limits temporal drifts. 

Finally, due to its limb-viewing geometry, the GPS RO has a higher vertical resolution (-1.5 km) 
than passive nadir sounders. This resolution is more comparable to radiosondes. Furthermore, the ratio 
between vertical and horizontal resolution (about 300 km) is consistent with that of quasi-geostrophic 
flows \Lindzen and Fox-Rabtnoi itz, 1989; hurstnski et al., 1997]. 


3.2. Configuration of a GPS Occultation 

Figure 1 shows the configuration of a GPS occultation. Three transmitters are used for positioning. 
One transmitter is used for clock correction. A fifth transmitter is used for the occultation. The rays 
from the source (occulted transmitter) traverse the Earth’s atmospheric limb to reach the receiver. 


Figure 1 




Figure 1. Radio occultation technique applied to GPS 


3.3. Bending Effect and Vertical Scanning 

Whenever an electromagnetic signal passes through the atmosphere, it is refracted according to 
Snell’s law, due to the vertical gradient of refractivity. This bending lengthens the ray path as compared 
with the straight path that would be obtained in a vacuum. Because the refractive index is not unity, 
the waves also travel slower than they would in a vacuum. The latter is more important m terms of 
delay than the bending. Finally, the GPS signal is affected by atmospheric absorption and scattering. 
For L-band (large) wavelengths, the scattering effect may be neglected with realistic liquid water and 
ice contents [Kursinski et ai, 1995]- 

At radio frequencies, it is not possible to make direct and accurate measurements of the refracted 
angle. However, if the two satellites are in relative motion, the refraction introduces a change in the 



Doppler shift of the received signal. After removal of the phase change due to the relative motion of the 
LEO with respect to the GPS, proper calibration of receiver and transmitter clocks, the extra phase 
change induced by the atmosphere can be isolated. 




Figure 3. Vertical scanning of the atmosphere during an occultation 

The overall effect of the atmosphere can be summed up by a total bending angle e and an 
asymptotic ray-miss distance p (commonly called impact parameter), as shown in figure 2. The distance 


figure 2 
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between the perigee point and the center of curvature is r^. The vertical scanning of the atmosphere 
is provided by the relative motion between the two orbiting satellites, as depicted in figure 3. Time 
dependency of both £ and p can be derived from accurate measurements of Doppler-shifted frequency. 


figure 3 


3.4. Inversion of Bending Angles to Refractivity 

The bending induced by the atmosphere ranges approximately between (2x10 at 80 km 
altitude and I" nearby the surface [A'ursmsfce et ai, 1997], The horizontal path covered by this bending 
is about 200~300 km long and 1.5~3 km wide, centered on the tangent point (see figure 4). A typical 


figure 4 


occultation lasts between 1 and 2 minutes. 

The variation of refractive index n along a limb path in the Earth’s atmosphere is dominated by 
the vertical gradient. To the first order, the refractive index field is spherically symmetric. If p is the 
impact parameter for the tangent ray whose radius is r, the bending angle £ can be expressed by . 



where x = rn [e.g. Eyre, 1994]. An elegant and direct way to obtain n from £ and p is to use Abelian 
transformation assuming a local spherical symmetry, i.e. . 


n (i;) 



( 2 ) 


[e.g. Eyre, 1994]. 

Errors in computing n by this approach result from (1) local spherical asymmetry, (2) non-coplanar 
rays, (3) non- vertical scanning (because both satellites drift during an occultation [Eyre. 1994]) and 
(4) an inaccurate upper boundary used to initiate the integral [Steiner et ai, 1999]. The integral 
formulation of the Abel transform spreads the errors in this boundary condition along the vertical 
and smoothes the observed bending angles variations. Consequently, the computed refractivity cannot 
reproduce the full sharpness of derived bending angles. 
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Figure 4. Bending of the rays during an occultation 


3.5. Interpretation of Refractivity 

For microwave frequencies, the refractivity N (defined by iV = 10® (n — 1)) of the atmosphere can 
be expressed by : 


N = 6i 


P . Pw I , L II/ 

y + * 2 :^ + f>3j2 + 


( 3 ) 


[e.g. Kursinskt et at.. 1997] where / is the frequency of the signal emitted by the transmitter in Hz, P 
the pressure of air (dry air and water vapor) in hPa, T the temperature in K, Pw the partial pressure in 
water vapor in hPa. Ne the electrons density in m“^, and W the particulates density (primarily liquid 
water) in g-m"®. The 6 , are constants; 6 i = 77.6 N-unit-K-hPa 62 = 3.73x10® N-unit-K hPa , 63 

= 40.3x10® N-unit-Hz- m®, 64 = 1.4 N-unit g"^ m®. 

The four refractivity terms in equation 3 are often refered to as (1) hydrostatic (or dry), (2) moist, 
(3) ionospheric and (4) scattering terms, respectively. The first term is due to the polarizability of 
atmospheric molecules. It is dominant below 60-90 km. The second term is due to the fact that the 
water molecule is polar. The ionospheric term results from the free electrons of the ionosphere. It can 
be removed by the use of the two GPS frequencies /i and [Kurstnski et ai, 1997]. The last term can 
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be neglected, as explained in section 3.3. Therefore, after removing the ionospheric component, iV can 
be approximated by [e.g. Bean and Dutton, 1968] ; 


\r L ^ I k 

^ = bl — + 02^2 ■ 


(4) 


4. Use of GPS for Atmospheric Retrievals and Data 
Assimilation 

There are commonly three levels of GPS radio occultation data available for data assimilation 
[Zou et ai, 1999]: (1) bending angles, (2) refractivity and (3) retrieved profiles of temperature or water 
vapor. To be assimilated, each level requires an estimate of the errors and an observation operator- 


4.1. Assimilation of Bending Angles 

These are considered as the “rawest” data (even though Doppler shifts themselves are even rawer 
data). An estimate of the error in each refracted angle measurement is needed. The advantage of using 
these data is that, at each point, the observation errors are mostly uncorrelated with respect to each 
other. 

However, the observation operator is quite complex [Eyre, 1994]. First, it has to interpolate 
the background information (typically gridded fields from a numerical weather prediction model) for 
temperature, humidity and pressure to each measurement location. This requirement stands for all 
approaches. Second, it has to compute the refracted angle 5 and its tangent linear model (or adjoint) for 
a given value of an impact parameter p. This implies computing the derivative of the local refractivity 
at each point, which is subject to high variability. 

Full ray-tracing codes have been developed by Mortensen and H0eg [1998] and Zou et al. [1999]. 
They both have the capability to reconstruct the ray path and simulate the overall bending angle from 
a given state of atmosphere. However, assimilating 30 profiles with a ray-tracing model in the National 
Centers for Environmental Prediction Data Assimilation System (NCEP-DAS) requires approximately 
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4 hours of CPU (Central Processing Unit) on a Cray C-90 computer [Matsumura et ai, 1999], 

Another approach consists of assuming spherical symmetry. Such one-dimensional operators have 
been developed by Healy [1998] and Palmer [1998]. The cost is significantly lower than a full ray-tracing. 
Like the 3D approach, it still involves local derivatives of refractivity. 

4.2. Assimilation of Retrieved Profiles of Refractivity 

In this approach, the bending angles have to first be inverted with an Abel transform to obtain 
profiles of refractivity as a function of altitude. The observation operator calculates refractivity from 
interpolated model variables using (4). This is much simpler and less computationally burdensome than 
simulating bending angles. 

The drawback to this approach is that additional errors are introduced as described in section 3.4. 
The main point is that errors are correlated along the vertical. These errors need to be estimated for 
accurate assimilation. 

4.3. Assimilation of Retrieved Profiles of Temperature and/or Humidity 

The retrieval error can be separated into three components [Rodgers, 1990] : (1) random error 
due to measurement noise (2) systematic error due to uncertain model parameters and inversion model 
bias (3) null-space error due to the inherent finite vertical resolution of the observing system (this is 
actually the background error component). The main point is that temperature and humidity errors in 
the retrievals are correlated both vertically and horizontally. These errors are difficult to estimate. 

The simplicity of such an approach is that the data to be assimilated are essentially the same as 
those used in a NWP model. Therefore, the observation operator is trivial. 

4.4. Summary on the Various Approaches 

Bending angles have the simplest error but the most complicated and expensive observation 
operator. Conversely, retrieved temperature or humidity profiles have the simplest observation operator, 
but the most complicated errors. The errors accumulate all the approximations made during the 
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conversion of raw GPS data into meteorological values. Assimilating refractivity represents a good 

alternative. It has a relatively low computing cost. 

The approach chosen here is to retrieve both temperature and humidity in a one-dimensional (ID, 
over the vertical) ‘off-line’ analysis based on a model forecast first guess and refractivity measurements. 
In the next chapter, the one-dimensional variational (IDV’.^R) analysis scheme is described. 

5. GPS-IDVAR Implementation 

5.1. Variational Theory 

The application of variational analysis to the retrieval of geophysical parameters has been discussed 
extensively by several authors [e.g. Rodgers, 1976, Eyre. 1993]. We try to minimize a cost function J 
with respect to a variable state of atmosphere x (state vector). This function is 

J[x] = {h[x) - y^VR-\h{x) - y°) + (X - x'^fR-^x - x"), (5) 

[e.g. Jazwinskt, 1970] where y° is the observation vector, h is the observation operator (non-linear), x 
is the background information (for example a 6-hour forecast). R is the error-covariance matrix of the 
observations, and B is the error-covariance matrix of the background. Hence h(x) is an estimate of the 
observations that would be made with a state of the atmosphere x. We assume the background and 
observation errors to be unbiased and uncorrelated with each other. 

The minimum variance problem can be solved using a quasi-Newton iteration, re., 

Xi+i x^ + {H'[R-^Hi + HJ R-^ [y° -h[x) + Hi (x; - x*-)) , (6) 

[e.g. Rodgers, 1976] where the subscript i denotes the iteration number. Hi is the tangent linear model 
of the observation operator h. or Jacobian, and Hj the adjoint. We define convergence as the iteration 
at which the quantity J [x.] has not changed by more than 2% from the previous iteration. In the 
dataset we used, this distance at convergence is in the mean equal to about two times the number of 
GPS observations used to resolve one occultation. Healy et al. [2000] and Palmer et al. [2000] have 
found this parameter can be used to quality control the GPS observations. 
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The state vector contains conventional meteorological variables on several pressure levels. After 
convergence, it represents a solution (the analysts) that has an optimal position with respect to both 
the observations and the background information. These distances, formally expressed in J, are 
normalized and weighted by the inverse of the matrices R and B, respectively. They represent physical 
constraints that help choose a solution among the infinite number of possible atmospheric temperature 
and humidity profiles that would match the observations. 

5.2. The State Vector 

The state vector includes the temperature and negative of the log of the specific humidity at the 
background pressure levels. The sea level pressure was added to the state vector as a supplementary 
degree of freedom. To reduce computational cost, the amount of humidity is taken as equal to a very 
small value (10"^) and removed from the state vector above a specified level (50 hPa). 

In the implementation, the state vector extends only to the lowest perigee point of the occultation. 
If n is the total number of levels, and np„ the number of model levels with water vapor m the state 
variable, the state vector x contains (nr = n + np„, + 1) elements : 

X = (Tl , ^2 - log - log 92 log 9np„ , P,eale«el ) ■ (0 

5.3. Background Information 

We used two sets of background information from the Goddard Earth Observing System (GEOS) 
model for our analysis. Both included upper stratospheric levels (GEOS-Strat). One set was a 6-hour 
forecast. Available archived data were specified on 18 selected pressure levels (surface up to 0.4 hPa). 
The other set was the ’assimilation' product on the 46 a levels of the General Circulation Model (GCM) 
from the surface up to 0.1 hPa. The Data Assimilation Office Data Assimilation System (DAO-DAS) 
uses an Incremental Analysis Update (lAU) approach [Bloom et al, 1996]. It consists of applying 
gradual increments within the 6-hour assimilation window. The result of this operation in the middle of 
the window is called the ’assimilation’. So the assimilation background we used already included some 
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information from the observations (especially radiosondes). 

The background error-covariance matrix uses variances from Joiner and Rokke [2000]. For 
temperature, it includes inter-level correlations through the following formula : 


Cov{Ti,Tj) = exp"^ ^ ^ exp 


( 8 ) 


where Cov{Ti,Tj) expresses the temperature error covariance between two levels i and j, (Tt. is the 
temperature standard deviation for the pressure level P,, T,- is the temperature at that level, AL = 0.1 
and AT = 30 are constants which were empirically adjusted. Two levels close in log of pressure and 
temperature will be considered as having highly-correlated errors. Also, in order to take into account 
the variable altitude of the tropopause, we added a feature to the scheme to fit the local maximum 
of model error to the altitude of the minimum of temperature. Figure 5 shows the diagonal of the 
background error covariance matrix for temperature. 

The assumed errors are relatively small. However, we did an experiment in which we multiplied 
the temperature errors by a factor 1.4. The results obtained by comparison with radiosondes (as shown 
later) were unchanged. The cost function J at convergence was increased. 

The values for humidity variances are listed in Table 1. They were slightly modified from Joinei 
and Rokke. [2000] after optimization by experiments. For the inter-level correlation only the term 
representing the exponential decrease related to the vertical distance between two levels in (8) is taken 
into account. 

We assume a sea level pressure standard deviation of 2.5 hPa. Background temperature, humidity, 
and sea level pressure errors are assumed to be uncorrelated. 


Figure 5 


Table 1 
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ADAPTATION OF THE B MATRIX TO THE TROPOPAUSE 


Temperature in K (dash line) 

180 200 220 240 260 280 300 320 



Background Standard Deviation (K) 


Figure 5. Assumed temperature errors for the background : plain line : for a tropopause at 250 hPa 

(typically mid and high latitudes) - diamonds : for the temperature profile with dashes (latitude 7°N). 
Table 1. Assumed log of specific humidity errors for the background. 


Pressure (hPa) 

1000 

850 

700 

500 

400 

300 

250 

200 

150 and beyond 

Log(q) std. dev. 

0.40 

0.40 

0.40 

0.20 

0.25 

0.30 

0.45 

0.60 

0.70 
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5.4. The GPS Observations 


During the GPS/MET mission, there were some time periods when the Anti-Spoofing (A/S) 
encryption was turned off by the USDoD, We used the data from one of these, called ’Prime Time 2’ 
from June 21st to July 4th 1995. The dataset, a total of 797 occultations, is shown in figure 6. 


figure 6 


We assumed the following errors in the refractivity data : 1% below 5 km, 0.2% up to 30 km. 
They follow the estimations of Kursinski et al. [1997]. Above 30 km. several sources of error, negligible 
in the troposphere, become gradually more important as the refractivity becomes very small [Kursmskt 
et al., 1997]. Hocke, [1997] identified “wavelike structures” in the upper stratosphere in GPS retrieved 
temperatures. Since we cannot determine whether these are real or simply noise in the observations, we 
simply chose not to trust the GPS refractivity observations above 30 km. We therefore assigned a 50% 


error in observed refractivity above that altitude. 


GPS/MET 

June 21 - July 04 1995 
Locotions of occultations 

+ OOH A 06H 
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Figure 6. GPS/MET Occultations with no GPS encryption, June-July 1995, processed by JPL. The 
chart is one map per day and the dates are shown at the bottom of each map, with the number of 

occultations for the day in parentheses. 

5.5. The Observation Operator 

The observation operator converts temperature, humidity and pressure profiles into refractivity 
profiles e.xpressed as a function of altitude. It contains the physics of the measurement and an 
appropriate space-time interpolation. We assume here that all observations happen at the synoptic 
time for which our background estimate is valid. The background profile at the observation location 
is obtained by interpolating bilinearly between the model grid points. It follows that (1) computing 
refractivity and (2) mapping pressure levels onto altimetric levels are the two main features of this 

operator. 

Instead of performing the analysis on the levels of the forecast, it is possible to work on more 
levels, for example the levels of the observations. It seems more appropriate to do so in order to obtain 
analyzed profiles that account for the high resolution features detected by the GPS and that cannot be 
tracked with a coarse analysis. However, the profiles are to be assimilated in a global circulation model 
(GCM) and thus must be brought back to a lower vertical resolution, more consistent with the model s 
dynamics. This is the reason why we chose to work on the levels of our background information. 
Moreover, the cost of computing on 18-46 levels is much smaller than 60-90 (typical number of GPS 
refractivity observations for one occultation). 

There are two different ways of performing the forward operation : (1) computing refractivity 
at each model level, then extrapolating refractivity profile expressed in pressure levels into a profile 
expressed in altimetric levels ; (2) extrapolating model values from each pressure level to the altitudes 
of the GPS observations, then computing at each point the refractivity value. 

Since our goal was to obtain the best profiles for further assimilation, we decided to calculate 
refractivity first and to interpolate it after. This requires only one interpolation. In addition, calculating 
refractivity with interpolated values for temperature and water vapor content is likely to generate 
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complexity in the computation of the Jacobian. 

Refractivity values are computed for each model level. For each of these levels, the allimetric 
altitudes are calculated using a hydrostatic integration. Assuming Pyeaievei corresponds to 0 meter 
altitude, the refractivity values are interpolated to the altitudes of the observations using a linear in log 
of pressure interpolation scheme. 

5.6. Linearized Version of the Observation Operator 

In (6) we must not only have an observation operator /i, but also it s linearized version H . H is 
the Jacobian or partial derivatives of the observation operator with respect to the various elements 
of the state vector. We derived an analytical formulation for // . As a check, it was compared with a 
computation by finite differences as in [Eyre et ai, 1993]. 

The refractivity at one level is found to be sensitive to a variation in temperature and humidity 
at that level and under that level. This is the result of the integrated altitude calculation using 
hydrostatic equilibrium. A change in the state vector at a given altitude has no influence on the 
simulated refractivities located below, but influences the altitudes of the pressure levels above, and 
hence the refractivity. Figure 7 shows a few columns of the H matrix. 

An increase in temperature means a decrease in refractivity at the same level. It has no influence 
on the calculated refractivities below. The value of the refractivity at each pressure level located above 
does not change. However, due to the increase in temperature and hydrostatic integration, all the 
pressure levels above are put higher. This results in shifting upwards the upper part of the refractivity 
profile, and increasing refractivity for a given altitude. 

Increasing humidity increases the local refractivity (the plot is for ~log(q)). It has no significant 
influence through the hydrostatic integration. 

The sea level pressure has no direct influence on the refractivity values of each model levels. 
However, it increases the pressure difference between each model level and the sea level, which is the 
same as increasing the altitudes. The final effect is hence to increase refractivity for a given altitude. 


Figure 7 
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-1.0 -0.5 0.0 0.5 1.0 

Sensitivity of refroctivity to temperature, humidity 
and sea level pressure (normalized) (1995 — 06—29 02:1 8Z) 

Figure 7. Jacobian for the GPS IDVAR. Curves for temperature at 607 and 94 hPa, for humidity at 
303 hPa, and sea level pressure. 

6. Theoretical improvement of GPS data over background 
information and Monte-Carlo simulation 


Knowing the errors contained in the GPS measurements, it is possible to infer the theoretical 
improvement brought by these data even before having real GPS observations. Rodgers [1990] has 
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proposed linear methods to evaluate the accuracy of a retrieval. They do not require simulated 
observations, contrary to non-linear methods. 

6.1. Linear Analysis 

Besides solving the temperature/humidity “ambiguity problem” the IDVAR approach gives also 
an appraisal of the quality of the analysis. Under the hypothesis of linearity and Gaussian error 
distributions, the error-covariance matrix of the analysis is 

P“ = + B~^y\ (9) 

[e.g. Rodgers, 1990]. The elements on the diagonal represent the variance of the solution with respect 
to the True state’ that we try to estimate. The interpretation of the results given by this method must 
be taken with caution, because of the hypotheses used: e.g. no bias, no correlation between different 
types of variables, perfectly estimated background and observation errors. 

In the calculations, we assume the same errors in the background for all latitudes. We assume a 
GPS refractivity vector containing 56 elements, from 1 to 56 km altitude. The associated errors are 2% 
below 5 km and 0.2% above (up to the top). To gauge the influence of the latitude, we use the same 
atmospheric profiles as Joiner and da Silva [1998]. 

Figure 8 shows the estimated standard deviation at various levels for both the background and the 
analysis (see legend as “Linear analysis”), for an atmospheric profile located at 30°. At lower latitude 
(18°, see figure 9), below 400 hPa, the water vapor content is important enough to play a role via the 
wet term in the refractivity in (4). But at higher latitudes (37° and 63°. see figure 10), the influence of 
water vapor is less important. The retrieval is therefore not as accurate in deriving fractionnal water 
vapor content. Concerning the temperature, the analysis has a significantly lower standard deviation 
than the guess. This result stands for all latitudes. 


Figure 8 


figure 9 
figure 10 
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6.2. Monte-Carlo Simulation 

In order to confirm the previous results and to test the IDVAR, a fully non-linear Monte-Carlo 
study is performed. Four experiments were carried out. 

In all experiments, 1000 atmospheric profiles are created. A basic atmospheric profile (the “truth”) 
is perturbed. The direction of the perturbations originate from the eigenvectors of the background 
covariance matrix. In each orthogonal direction, the intensity of the perturbation is obtained multiplying 
the corresponding square root of the eigenvalue by a random number (normal distribution). Finally, 
the sum of all the perturbations and the original “true state” give the background estimate. 

6.2.1. Observations : No Error Added. For the observations, the refractivity profile is 
taken as the exact result of the observation operator applied to the true state. No error is added. We 
ensure that 1000 profiles are enough to generate a distribution with no significant bias as assumed in 
the theory (e.g. after generation of the population the temperature bias is less than 0.06K). Figure 8 
shows the statistical results of the experiment. They are better than what is expected from the linear 
error analysis which assumed observational errors. We can see that by neglecting the observation errors 
in an OSSE, we may slightly overestimate the impact of the observations. 

6.2.2. Observations : Errors Added. Here, the observed refractivity profiles are perturbed, 
are added about the exact simulated refractivity obtained from the true state, using the 

eigenvalues and eigenvectors of the observation error covariance R. Three experiments are made 
for three different atmospheric profiles, each of them characteristic of a specific latitude (Tropics. 
Mid-latitude, High-latitude). 

Figure 10 and figure 1 1 show a good agreement with the linear error analysis results. This yields 
an important point : the IDVAR works correctly with simulated data (total of 3000 analyses, all 
converged). 

Under the above hypotheses, refractivity data are expected to bring some improvement for 
temperature at all latitudes. For the humidity, the best improvements may happen in regions with a 
relatively high amount of water vapor (i.e. Tropics, mid to lower troposphere). Palmer et al. [2000] 


figure 11 
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and Healy et al. [2000] have obtained similar results. 




Figure 8. Standard deviation of : top : temperature ; bottom : log specific humidity. Linear analysis 
assumes noise in the GPS observations; results are in plain line for the guess estimate and short dashes 
for the analysis. Monte-Carlo simulation is performed with no noise in the observations: the curves are 
dotted line for the guess and long dashes for the analysis. 
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pressure 

Figure 11. Same as figure 9, latitude 18®, but for the temperature. 

6.3. Averaging Kernels 

The repartition of information brought by the GPS IDVAR can be made by looking at the 
averaging kernels [Rodgers, 1990]. They are computed the following way: (1) for the temperature, a 1 K 
impulse is added in the state vector at one level (2) the observation operator is applied to that state of 
atmosphere to obtain a refractivity profile which contains effectively the 1 K impulse (3) we perform a 
IDVAR analysis using for the background the unperturbed state vector (4) the difference between the 
retrieved temperature and the original state vector gives the response of the system. This is repeated 
for each level of the guess estimate. If the observations contained perfect information content and no 
weight was given to the background, the response would be equal to the original unit impulse, and 
focused at the same level only. Figure 12 testifies that the IDVAR analysis keeps the relatively high 


Figure 12 
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vertical resolution of the observations. 

In the given example, the observations start around 5 km altitude. The state vector is shortened 
so that only levels located above that altitude are modified by the analysis. Consequently, the response 
of the system to a IK impulse located anywhere below that level is zero. 

GPS 1DVAR AVERAGING KERNELS 

35 km 
30 km 

25 km 

-g 20 km 

D 

^ 1 5 km 

1 0 km 

5 km 
0 km 

-0.2 0.0 0.2 0.4 0.6 0.8 

in Kelvin/Kelvin 



Figure 12. Temperature averaging kernels for the GPS IDVAR. 


7. Sensitivity of the retrievals to various parameters 

7.1. Sensitivity to the Number of Levels in the Guess 

We will now focus on a single occultation which occurred on June 21st 1995, at 00:03UTC, at 
(10®S, 168^ W). Three relevant temperature profiles are discussed here : the first-guess estimate (6-hour 
foreccLst of the GEOS model), the analysis obtained through the IDVAR retrieval, and the direct 
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retrieval (assuming a humidity profile). The latter depicts how the “raw” refractivity data see the 
atmosphere. Here we assume that the temperature/humidity ambiguity problem will not lead to a 
misinterpretation of the data. The refractivity data start at 6600 m for this observation. Therefore, 
the risk of misinterpretation is limited, even though the occultation occurred in the Tropics. For the 
IDVAR analysis, the error in refractivity is taken as 2.% up to 5 km, and 0.2% above, with no inter-le\el 


correlation. 

Occultolion: Lot -10.21 deg. Lon 


-167.6deg. 1995 06 21. 00:03UTC, Tronsmitter 28 



Temperolure (K) 


Figure 13. Temperature profile, IDVAR analysis performed on 18 levels. 

7.1.1. Analysis on 18 levels. Figure 13 shows the result of an analysis performed with 



background information on 18 levels. These levels are the available archived levels. The IDVAR works 
on the same levels as those of the background and hence the retrievals are also on 18 levels. 


The fine resolution of the GPS measurements as inferred by the direct retrieval appears to produce 


some high frequency structure in the stratosphere. The IDVAR tries to simulate this structure, as the 
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observation errors are relatively small as compared with the projected background errors. However, 
due to the coarse resolution of the state vector (bcised on the 18 levels of the guess), the analysis cannot 
fully resolve these structures. 

The background used to perform the error analysis study and Monte-Carlo simulation was on 
18 levels. These studies did not take into consideration any inconsistency between the resolution of the 
observing system and the guess. However, the data have a finer resolution so that a coarse- resolution 
IDVAR cannot fully resolve the structure. 

7.1.2. Analysis on 46 levels. Another first-guess estimate is used here. This guess consists of 
46 cr levels, with a top at 0.1 hPa. The a levels are converted to pressure. Figure 14 shows the result of 
the analysis of the same refractivity profile on 46 levels. We can see some “waves'’ in the analysis at 
this resolution. We also see a large increment at the top of the profile. Note that the .structure at the 
tropopause is now much better resolved than before at low resolution. This shows the need to analyze 
the GPS at a higher resolution than the GEOS-Strat analysis if the full impact is to be realized. As 
explained in section 5.4, we will subsequently increase the errors of GPS observations in the analysis 
above 30 km. 


Figure 14 
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Occullatlon: Lot -10.21 deg. Lon -167.6deg. 1995 06 21. 00:03UTC. Tronsmilter 28 



Figure 14. Temperature profile, IDVAR analysis performed on 46 cr levels. 

7.2, Sensitivity to the gravitational force 

One unexpected result is the sensitivity of the analysis to the gravitational constant used in the 
analysis. We used two formulations for g. (1) ^=9.80665 m*s““ (2) g=g[z) [Healy, 1998]. The second 
gave better results by comparison with radiosondes. 

Depending on the approximation used for the definition of g, differences of order IK can appear. 
This sensitivity to the gravitational constant is specific to GPS measurements. It arises from the 
conversion from pressure (the levels of the analysis) to altitude (the levels of the observations). 
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7.3. Sensitivity to the Surface Pressure 

The surface pressure also affects the computation of the altitudes. It is used as a reference to 
integrate the hydrostatic relationship. 


The analysis increment is defined as the difference between the IDVAR retrieval and the initial 
background. Figure 15 shows two different increments. The first one is calculated with a fixed reference 
pressure, i.e. the same as in the background estimate. In the second implementation, the reference 
pressure becomes part of the state vector and hence can be updated. Overall, radiosonde comparison 
statistics are improved by this modification. When the reference pressure can be updated, the whole 
refractivity profile can be moved upwards or downwards. When it is not possible for the system to do 
that, its only alternative is to increase or decrease the temperature on several layers so that the layers 
can get thicker or thinner. 

Adding this degree of freedom enables the IDVAR to account for any inconsistencies that may 
exist between the guess and the observations concerning the definition of the sea level. It may also limit 
the errors due to the interpolation of the sea level pressure from the background fields. 

To illustrate this, taking the example of occultation 1995-06-28 at 09:31Z, the computed altitude of 
the level lOhPa is moved from 31,160 to 31.290 meters (after minimization) when the sea level pressure 
(SLP) cannot be moved by the analysis (i.e. the difference is created by the adjustment made on the 
temperatures). When the SLP is included in the state vector, the temperature increments are smaller, 
but the computed altitude is about the same as before at 31,280 meters. 


Figure 15 


8* Comparisons with Radiosondes 

8.1. Collocation 

The statistics shown here are computed by collocating results of GPS IDVAR analyses with nearby 
radiosondes (RS). The criteria for collocation are +/- 3 hours, less than 280 km. The temperature 
(humidity) direct retrievals obtained using an estimate of humidity (temperature) from the European 




OccuUoUoo: Lot 49.65d#9, Loo 24.710eg. 199S 06 28. 09:31UTC, Tron»miU«r 18 



Anolysis increment in temperature (K) 

Occultotloo: Lot 49.65d*9, Loo 24.71 deg. 1995 06 28. 09:31UTC. TronimUter 18 



Analysis increment in temperature (K) 

Figure 15. Temperature analysis increment. Above : sea level pressure of the background assumed as 
correct - Below : sea level pressure included in the state vector 
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Centre for MediunvRange Weather Forecasts (ECWIWF) analyses are also shown [hursinski et aL, 

1997]. 

Among the 797 GPS occultations, about 150 match the collocation criteria. Results are not shown 
for the southern hemisphere (beyond 30^S) because of a very small number of comparisons. Only one 
IDVAR minimization process did not converge. As a background check, collocations with discrepancies 
more than 5 K between background and RS are removed from the statistics. It is necessary to apply 
some quality control to the RS and the direct retrievals. Therefore, we also remove discrepancies more 
than 5 K between the RS and the direct retrieval. No quality control is performed on the IDVAR 
analyses themselves. 

The same policy is applied for humidity, expressed in terms of difference of specific humidity (in 
percent). All differences, background minus RS or RS minus direct retrieval, greater than 100% are 
removed. 

8.2. Temperature 

Figure 16 shows a single temperature profile which exhibits the advantages of combining model Figure 1 
estimates with real observations. (1) Both direct retrieval and IDVAR analysis are able to resolve the 
tropopause better than the GEOS assimilation. The background has a warm bias at the tropopause 
of about 5K. (2) In the lower layers, the analysis follows the background, whereas the direct retrieval 


notably diverges from the RS. 
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Occultation: Lot 58,94deg, Lon -*65,88deg. 1995 07 03* 12:21UTC, Transmitter 16 



Figure 16. Temperature profile showing the advantages of using a variational approach with GPS data. 

Figure 17 presents the biases and root-mean square (rms) difference from collocated radiosondes 
for temperature in the northern hemisphere (above 30®N). The chart also indicates (to the right) the 
number of data used to calculate each point. That number decreztses towards the surface because most 
of the occultations did not reach the surface. The balloon bursts above the tropopause explain the 
decrezLsing number of comparisons in the stratosphere. 

A significant temperature bias present in the background data at the tropopause is removed by the 


Figure 17 
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IDVAR analysis, between 250 and 100 hPa. The standard deviations are significantly reduced above. 
Below 500 hPa, a slight degradation can be detected. Overall, the rms of the GPS IDVAR is reduced 
as compared with the background at most altitudes. 

Similar results are obtained in the Tropics. Figure 18 shows the temperature bias in the tropics 
(between 30*S and 30‘’N). Around 100 hPa, the GPS IDVAR reduces a bias even more significant 
than in the northern region. Figure 19 demonstrates a net reduction of the rms of about l.OK at the 
tropopause. 

Figure 20 represents the increments, i.e. the difference between IDVAR retrieved temperature and 


Figure 18 


Figure 19 


Figure 20 


the initial value given by the GEOS model. They are largest at the tropopause. 
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Figure 17. Northern hemisphere, comparisons with radiosondes temperatures (+/- 3 hours, less than 
280 km) for the GPS ID VAR, the background and the direct retrieval. Curves show bias and rms. 
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Figure 20. GPS IDVAR temperature increments. 


8.3. Humidity 

Many radiosondes we used either did not report the humidity, or the reported humidity was 
different from the GEOS background by more than 100%. It is known that humidity present smaller 
scale features than temperature. There are also known problems with RS humidity meeisurements 
[Soden et aL, 1996]. The collocation criteria chosen may be responsible for many of the discrepancies. 
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This results in a small number of comparisons. We were not able to find a sufficient number of 
radiosondes presenting a reasonable difference with the background in the tropics. For this reason, only 


the northern hemisphere results are presented. 

Figure 21 shows bias and rms with respect to collocated RS. Only the 400 hPa level shows an 
improvement over the background. Some degradation of the IDVAR as compared with the background 
is shown at other levels. 


Figure 21 
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Figure 21. Northern hemisphere, comparisons with radiosondes humidity. Curves show bias and rms. 
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8.4. Observation residuals 

Residuals are defined as the difference between the observation and the simulated refractivity 
calculated from the model background. Since refractivity spans values of several orders of magnitude, 
we are looking at it in terms of percent of observed refractivity. 

The refractivity varies locally primarily with the inverse of temperature. The fact that the 
background has a warm bias at the tropopause implies a local negative bias in refractivity as compared 
to the observations. However, since the above altitudes are computed using the temperature at the 
levels below, all the refractivity values above don’t move. However, the profile of geometric altitudes 
corresponding to the same pressures is shifted upwards. This results in an increasing positive refractivity 

bias (see figure 22). 

In the tropics (figure 23), the negative bias due to the missed tropopause is located higher, around 
1.5 km. Below 5 km, the refractivity differences grow significantly, indicating that (1) phenomena 
are not simulated correctly by our observation operator and/or (2) there are important errors in the 
measurements. Horizontal gradients and multipath may appear in either category, depending on the 
way the analysis is made (i.e. ID or 3D, refractivity or bending angle, geometrical optics or Fresnel 

inversion). 

In the IDVAR theory, we assumed that both observations and the background are not biased. 
Clearly they are biased. An interesting point is that even with a bias, the IDVAR removes some of that 
bias : the curve for the bias of the IDVAR refractivity residual is close to zero. However, if we look 
carefully, some bias remains in the residuals. This is most apparent near the tropopause in the tropics 
where the background bias was the strongest. 


figure 22 
figure 23 
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Percents of observed refroctivity value 
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Figure 22. Difference between observed refractivity and simulated refractivity (from the GEOS back- 
ground), in percents of the observed refractivity, in the northern hemisphere. 
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Figure 23. Same as figure 22 (observation residuals), but for the tropics. 


9. Conclusions and Future Directions 

An overview of the various techniques possible to assimilate GPS data has been presented. We 
implemented a IDVAR analysis of refractivity. This presents a good trade-off between computational 
cost and accuracy. The one-dimensional observation operator is efficient and accurate. Temperature 
improvements are both expected and observed with GPS/MET data. We saw improvement at some 
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altitudes in the analysis of humidity with GPS data, but this was accompanied with some degradation 
at other altitudes. 

There are assumptions made when deriving refractivity from bending angles (i.e. local spherical 
symmetry). Even with this assumption, we have shown that significant improvement in temperature 
can be obtained by combining GPS observations with an accurate background. However, the limited 
impact in regions with high water vapor amounts (i.e. below 500 hPa), the strong discrepancy between 
simulated and observed refractivities and the unexpected sensitivity of the retrievals to the formulation 
of the gravitational constant testify that some phenomena still need to be studied with care before 
making any statement related to global trends or impact on NWP based on GPS observations. 

After showing a positive impact on the analysis, the next step to demonstrate the usefulness of 
GPS in weather prediction is to run an impact experiment. The small amount of data collected daily 
with GPS/MET currently limits the significance of such an experiment, regardless of the quality of the 
data. More GPS data with improved receivers are expected in the future. They will give an opportunity 
for more global and case studies. 
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